//#include <iostream.h>
//#include <iomanip.h>

// ----------------------------------------------------------------
// Alexander Doxiadis 
// adox@nikhef.nl
// 13-09-2010 
//
// simple function that takes as input two histogram (same binning,
// same variable): 1st is the one with loose lepton cuts and 2nd
// with tight lepton cuts, both before b-tag
//  
// for MM details see:
// https://svnweb.cern.ch/trac/atlasgrp/browser/Physics/Top/TTBar2010/Note5/note5_July29.pdf
//
// the function needs the following flags:
// 
// IJet == 1: 1  jet , pt > 20
// IJet == 2: 2  jets, pt > 20
// IJet == 3: 3  jets, pt > 20
// IJet == 4: 4+ jets, pt > 20
//
// IJet == 5: 1  jet , pt > 25
// IJet == 6: 2  jets, pt > 25
// IJet == 7: 3  jets, pt > 25
// IJet == 8: 4+ jets, pt > 25
//
// IBtag == 0 : returns estimate without b-tag
// IBtag == 1 : returns estimate after b-tag
// 
// IError == 0 : statistical uncertainty
// IError == 1 : full systematical uncertainty
// IError == 2 : full systematical + statistical uncertainty
// ----------------------------------------------------------------

TH1D* PlotQCD(TH1D* _loose , TH1D* _tight, Double_t eff_fake, Double_t eff_fake_error, Int_t IJet = 4, Int_t IBtag = 0, Int_t IError = 0){

  // Define the efficiencies:

  // whole 2010 data without btagging
  //Double_t eff_fake = 0.311631;
  //Double_t eff_fake_error = 0.00993165;

  // whole 2010 data with btagging
  //Double_t eff_fake = 0.264659;
  //Double_t eff_fake_error = 0.0274307;

  // 2010 data (2.9 pb^-1)
  //Double_t eff_fake = 0.294985;
  //Double_t eff_fake_error = 0.0247685;

  Double_t eff_real = 0.98;
  Double_t eff_real_error = 0.02;

  vector<Double_t> btags = btag_eff(IJet);
  Double_t b_eff = 0.;
  Double_t b_err = 0.;
  
  if (IBtag == 0){ // no b-tag cut
    b_eff = 1.;
    b_err = 0.;
  }
  else if (IBtag == 1){ // with b-tag cut
    b_eff = btags[0];
    b_err = btags[1];
  }

  // Get the plots:
  TH1D * loose = _loose->Clone("loose");
  TH1D * tight = _tight->Clone("tight");

  // Define the new plot:
  TH1D * qcd = tight->Clone("qcd");
  qcd->Reset();

  for(int i = 0 ; i < qcd->GetNbinsX() + 1; i++){
    
    Double_t N_loose = loose->GetBinContent(i); // NLoose
    Double_t N_tight = tight->GetBinContent(i); // NTight
    Double_t N_loose_error = sqrt(N_loose);
    Double_t N_tight_error = sqrt(N_tight);
    
    if (N_loose > 0){
      vector<Double_t> estimate = qcd_estimate(N_loose, N_loose_error, N_tight, N_tight_error, eff_real, eff_real_error, eff_fake ,eff_fake_error );

      cout << estimate[0] << endl;

      qcd->SetBinContent(i,b_eff*estimate[0]);
      if (IError == 0 ) qcd->SetBinError(i,b_eff*estimate[1]);
      else if (IError == 1) {
	Double_t sys_err = sqrt( 
				pow(b_eff*estimate[2],2)+           // systematic from efficiencies
				pow(b_err*b_eff*estimate[0],2)+     // systematic from 30% btag uncertainty + btag mc uncertainty
				pow(0.16*b_eff*estimate[0],2)+      // from closure
				pow(0.02*b_eff*estimate[0],2)+      // eff real from mc
				pow(0.195*b_eff*estimate[0],2)      // all other systematic like contamination and control region
				);
	qcd->SetBinError(i,sys_err);
      }
      else if (IError == 2) {
	Double_t sys_err = sqrt( 
				pow(b_eff*estimate[3],2)+           // systematic from efficiencies + stat uncertainty
				pow(b_err*b_eff*estimate[0],2)+     // systematic from 30% btag uncertainty + btag mc uncertainty
				pow(0.16*b_eff*estimate[0],2)+      // from closure
				pow(0.02*b_eff*estimate[0],2)+      // eff real from mc
				pow(0.195*b_eff*estimate[0],2)     // all other systematic like contamination and control region
				);
	qcd->SetBinError(i,sys_err);
      }
    }
    else {
      qcd->SetBinContent(i,0);
      qcd->SetBinError(i,0);
    }
  }  
  
  // Return the new plot:
  return qcd;

}

// ------------------------------------------
// function to return qcd estimate given 
// N_loose, N_tight and eff_real and eff_fake
// ------------------------------------------
vector<Double_t> qcd_estimate(Double_t N_loose, Double_t N_loose_error, Double_t N_tight, Double_t N_tight_error, Double_t eff_real, Double_t eff_real_error, Double_t eff_fake, Double_t eff_fake_error ){
  
  vector<Double_t> _error;
  _error.resize(4);

  // Final uncertainty
 Double_t Nf       = 0.00;
 Double_t Nf_error = 0.00;

 Double_t N_delta         = N_loose - N_tight;
 Double_t N_delta_error   = TMath::Sqrt(N_delta);

 // ----------------------------
 // Compute N_fake and its error
 // ----------------------------

 //-- N_fake
 Nf = ( eff_fake / (eff_real-eff_fake) ) * (eff_real*N_loose-N_tight);

 //-- uncertainty
 Double_t Nf_error      = 0.; 
 Double_t Nf_error_stat = 0.; 
 Double_t Nf_error_syst = 0.; 

 Double_t df_dEf = 0.; 
 Double_t df_dEr = 0.; 
 Double_t df_dNl = 0.; 
 Double_t df_dNt = 0.;
 Double_t df_dNd = 0.;


 // Nf = ( eff_fake / (eff_real-eff_fake) ) * [ (eff_real-1)*N_tight + eff_real*N_delta) ];
 
 df_dEf  = eff_real * (eff_real*N_loose - N_tight) / TMath::Power((eff_real-eff_fake),2); 
 df_dEr  = eff_fake * (N_tight-eff_fake*(N_loose)) / TMath::Power((eff_real-eff_fake),2); 
 df_dNt  =  (eff_real-1.) * ( eff_fake / (eff_real-eff_fake) ); 
 df_dNd  =  (eff_real)    * ( eff_fake / (eff_real-eff_fake) ); 
 
 Nf_error_stat = TMath::Sqrt(  
			     TMath::Power(df_dNd,2)*TMath::Power(N_delta_error,2)   + 
			     TMath::Power(df_dNt,2)*TMath::Power(N_tight_error,2) 
			     );
 
 Nf_error_syst = TMath::Sqrt(  
			     TMath::Power(df_dEf,2)*TMath::Power(eff_fake_error,2)  + 
			     TMath::Power(df_dEr,2)*TMath::Power(eff_real_error,2)
			     );

 Nf_error      = TMath::Sqrt(  
			     TMath::Power(Nf_error_stat,2) + 
			     TMath::Power(Nf_error_syst,2)
			     );
 
 _error[0] =  Nf;
 _error[1] =  Nf_error_stat;
 _error[2] =  Nf_error_syst;
 _error[3] =  Nf_error;
 
  return _error;
}


// -----------------------------------------
// function to return b-tag efficiency from
// qcd mc, 108180 r1306q + 108181 r1306r
// -----------------------------------------
vector<Double_t> btag_eff(Int_t IJet = 1){


  vector<Double_t> _btag;
  _btag.resize(2);

  Double_t untag = 0.;
  Double_t tag   = 0.;
  
  // exclusive jet bins, except 4+ jet, jet_pt > 20 GeV
  if (IJet == 1) {
    untag   = 322.;
    tag     =  22.;
  }
  else if (IJet == 2) {
    untag  = 123.8;
    tag    =  14.3;
  }
  else if (IJet == 3) {
    untag  = 36.3;
    tag    =  7.1;
  }
  else if (IJet == 4) {
    untag   = 16.0;
    tag     =  8.1;
  }
  
  // exclusive jet bins, except 4+ jet, jet_pt > 25 GeV
  else if (IJet == 5) {
    untag =  286.0;
    tag   =  20.0;
  }  
  else if (IJet == 6) {
    untag = 63.5;
    tag   = 8.2;
  }  
  else if (IJet == 7) {
    untag  = 13.5;
    tag    =  2.1;
  }  
  else if (IJet == 8) {
    untag =  8.6;
    tag   =  5.1;
  }
  
  Double_t syst = 0.30; // take the value of 30% uncertainty (c-tag uncertainty) from b-tag group
  
  _btag[0] = tag/untag;    
  _btag[1] = sqrt(pow(syst, 2) + (1. - _btag[0])/tag);
    
  return _btag;
  
}
